The field quantity $\vec{\mathcal{A}}(\vec{r},t)$ which is produced by the source $\vec{\mathcal{J}}_i(\vec{r},t)$ are functions of spatial coordinates ($\vec{r}$) and time ($t$).  Assume that these two independent variables can be separated as follows.
\begin{align}
	\vec{\mathcal{A}}(\vec{r},t) &= \vec{A}(\vec{r})T(t)\label{eqn:sepeqnrt}\\
	\vec{\mathcal{J}}_i(\vec{r},t) &= \vec{J}_i(\vec{r})T(t)\label{eqn:sepeqnrtT}
\end{align}

\emph{Question}?  Here we assumed that the independent variables could be separated by a multiplication operator.  Why multiplication?  Why not some other operator?  How about addition or maybe convolution, or maybe some kind of an inner product? This question may be answered by Green's functions.  See Balanis Ch 14 ADVANCED ENGINEERING ELECTRODYNAMICS.

Substituting (\ref{eqn:sepeqnrt}) into (\ref{eqn:waveeqn0}) yields,
\begin{multline}
T(t)\nabla^2\vec{A}(\vec{r})+\vec{{J}}_i(\vec{r})T(t)=\\\mu\varepsilon\vec{A}(\vec{r})\frac{d^2T(t)}{d{t}^2}+\left(\mu\sigma_e+\varepsilon\sigma_m\right)\vec{A}(\vec{r})\frac{dT(t)}{dt}+\sigma_m\sigma_e\vec{A}(\vec{r})T(t)\label{eqn:waveeqnsep}
\end{multline}
Notice that $\nabla^2$ only operates on $\vec{A}(\vec{r})$ and $\frac{d}{dt}$ only operates on $T(t)$. By dot multiplying both sides of (\ref{eqn:waveeqnsep}) with $\frac{\hat{a}}{AT}$, where $\vec{A}=A\hat{a}$ yields,
\begin{align}
\frac{\hat{a}(\vec{r})\cdot\nabla^2\vec{A}(\vec{r})}{A(\vec{r})}+\frac{\hat{a}(\vec{r})\cdot\vec{{J}}_i(\vec{r})}{A(\vec{r})}&=
\frac{\mu\varepsilon}{T(t)}\frac{d^2T(t)}{d{t}^2}+\frac{\mu\sigma_e+\varepsilon\sigma_m}{T(t)}\frac{dT(t)}{dt}+\sigma_m\sigma_e\label{eqn:waveeqnsep1}
\end{align}
Notice that the left hand side of (\ref{eqn:waveeqnsep1}) only depends on the spatial coordinate $\vec{r}$ and the right hand side only on the temporal variable $t$.  The only way for (\ref{eqn:waveeqnsep1}) to be true is for both sides of the equation be equal to a constant. Therefore, (\ref{eqn:waveeqnsep1}) reduces to the following differential equations,
\begin{align}
\frac{\hat{a}(\vec{r})\cdot\nabla^2\vec{A}(\vec{r})}{A(\vec{r})}+\frac{\hat{a}(\vec{r})\cdot\vec{{J}}_i(\vec{r})}{A(\vec{r})}&=\gamma^2=-k^2\label{eqn:waveeqnsepr}
\end{align}
\begin{align}
\mu\varepsilon\frac{T''(t)}{T(t)}+\left(\mu\sigma_e+\varepsilon\sigma_m\right)\frac{T'(t)}{T(t)}+\sigma_m\sigma_e&=\gamma^2=-k^2\label{eqn:waveeqnsept}
\end{align}
where, $(\frac{d}{dt}=')$ and $(\gamma^2=-k^2)$ is a constant. Equation (\ref{eqn:waveeqnsepr}) has the following solution,
\begin{align}
	\nabla^2\vec{A}(\vec{r})+\vec{J}_i=\gamma^2\vec{A}(\vec{r})=-k^2\vec{A}(\vec{r})\label{eqn:waveeqnsoln}
\end{align}
This can be verified by substituting (\ref{eqn:waveeqnsoln}) back into (\ref{eqn:waveeqnsepr}).
